In [1]:
include("program.jl")
using Plotly
importall DoubleIntegralUtils
importall DoubleIntegralNewtonCotes
importall DoubleIntegralGauss
importall DoubleIntegralPadua

Plotly javascript loaded.

To load again call 

init_notebook(true)

In [2]:
function draw_functions(a,b,c,d,precision, function_list)
    function get_scatter(f,a,b,c,d,precision)
        domain3d = [(x,y) for x in linspace(a,b,precision) for y in linspace(c,d,precision)]
        return Plotly.scatter3d(x=[t[1] for t in domain3d],y=[t[2] for t in domain3d],z=[(f(t[1], t[2])) for t in domain3d], mode="lines")
    end
    
    Plotly.plot([get_scatter(func,a,b,c,d,precision) for func in function_list])
            end

WARNING: Method definition (::Type{PlotlyJS.SyncPlot{TD} where TD<:PlotlyJS.AbstractPlotlyDisplay})(PlotlyJS.Plot{TT} where TT<:PlotlyJS.AbstractTrace) in module PlotlyJS at /Users/wektor/.julia/v0.6/PlotlyJS/src/display.jl:136 overwritten at /Users/wektor/.julia/v0.6/PlotlyJS/src/displays/ijulia.jl:187.
Out[2]:
draw_functions (generic function with 1 method)
In [3]:
function test_function(f,a,b,c,d,exact_value,n)
    ns = [2:1:n;]
    Plotly.plot([Plotly.scatter(x=ns, y=[abs(composite_double_trapezoid(f,a,b,c,d,i,i)-exact_value)/abs(exact_value) for i in ns], name="Composite Trapezoid Method Error"),
                 Plotly.scatter(x=ns[1:2:end], y=[abs(composite_double_simpson(f,a,b,c,d,i,i)-exact_value)/abs(exact_value) for i in ns[1:2:end] ], name="Composite Simpson Method Error"),
                 Plotly.scatter(x=ns[1:1:min(n-1,14)], y=[abs(double_n_by_m_order(i,i,f,a,b,c,d)-exact_value)/abs(exact_value) for i in ns[1:1:min(n-1,14)]], name="Newton-Cotes Cubature Error"), 
                 Plotly.scatter(x=ns, y=[abs(double_gauss_chebyshev((x,y) -> (move_function(f,a,b,c,d)(x,y))*sqrt(1-x^2)*sqrt(1-y^2),i,i)-exact_value)/abs(exact_value) for i in ns], name="Gauss Cubature Error"),
                 Plotly.scatter(x=ns, y=[abs(padua(move_function(f,a,b,c,d), i)-exact_value)/abs(exact_value) for i in ns], name="Padua Cubature Error")])
end

function test_exact_digits(f,a,b,c,d,exact_value,n)
    ns = [2:1:n;]
    Plotly.plot([Plotly.scatter(x=ns, y=[exact_digits(composite_double_trapezoid(f,a,b,c,d,i,i),exact_value) for i in ns], name="Composite Trapezoid Method Exact Digits"),
                 Plotly.scatter(x=ns[1:2:end], y=[exact_digits(composite_double_simpson(f,a,b,c,d,i,i),exact_value) for i in ns[1:2:end] ], name="Composite Simpson Method Exact Digits"),
                 Plotly.scatter(x=ns[1:1:min(n-1,14)], y=[exact_digits(double_n_by_m_order(i,i,f,a,b,c,d),exact_value) for i in ns[1:1:min(n-1,14)]], name="Newton-Cotes Cubature Exact Digits"), 
                 Plotly.scatter(x=ns, y=[exact_digits(double_gauss_chebyshev((x,y) -> (move_function(f,a,b,c,d)(x,y))*sqrt(1-x^2)*sqrt(1-y^2),i,i),exact_value) for i in ns], name="Gauss Cubature Exact Digits"),
                 Plotly.scatter(x=ns, y=[exact_digits(padua(move_function(f,a,b,c,d), i),exact_value) for i in ns], name="Padua Cubature Exact Digits")])
end
Out[3]:
test_exact_digits (generic function with 1 method)

Ze względu na limity dotyczące rozmiaru plików, zachęcam do wykorzystania zakomentowanego kodu.

Przykłady wykorzystania programu.

Funkcja move_function(f,a,b,c,d) zwraca taką funkcję g, że $ \int_a^b \int_c^d f(x,y) dydx = \int_{-1}^{1} \int_{-1}^{1} g(x,y) dydx $.

In [4]:
f(x,y) = 1/(1+x^2+y^2)
g = move_function(f,0,1,0,1)
gf = (x,y) -> g(x,y)*sqrt(1-x^2)*sqrt(1-y^2)
Out[4]:
(::#38) (generic function with 1 method)
In [5]:
@show composite_double_trapezoid(f,0,1,0,1,10,10)
@show composite_double_simpson(f,0,1,0,1,10,10)
@show Float64(double_n_by_m_order(10,10,f,0,1,0,1))
@show double_gauss_chebyshev(gf, 10,10)
@show padua(g, 10)

composite_double_trapezoid(f, 0, 1, 0, 1, 10, 10) = 0.6388700512475036
composite_double_simpson(f, 0, 1, 0, 1, 10, 10) = 0.6395105575477145
Float64(double_n_by_m_order(10, 10, f, 0, 1, 0, 1)) = 0.6395103671022478
double_gauss_chebyshev(gf, 10, 10) = 0.6436706081080388
padua(g, 10) = 0.6395103516476914
Out[5]:
0.6395103516476914

Testy funkcji $\frac{1}{1+x^2 + y^2}$.

In [6]:
f(x,y) = 1/(x^2 + y^2 + 1)
draw_functions(0,1,0,1,100, [f])
Out[6]:
In [7]:
exact_int = 0.639510351870311001962693085427323679
Out[7]:
0.639510351870311
In [8]:
test_function(f,0,1,0,1,exact_int, 30)
Out[8]:
In [9]:
test_exact_digits(f,0,1,0,1,exact_int, 30)
Out[9]:

Testy funkcji Franke'go.

In [10]:
f1(x,y) = 0.75 * exp(-(9*x-2)^2/4 - (9*y-2)^2/4)
f2(x,y) = 0.75 * exp(-(9*x+1)^2/49 - (9*y+1)/10)
f3(x,y) = 0.50 * exp(-(9*x-7)^2/4 - (9*y-3)^2/4)
f4(x,y) = -0.2 * exp(-(9*x-4)^2 - (9*y-7)^2)
f(x,y) = f1(x,y) + f2(x,y) + f3(x,y) + f4(x,y)
exact_int = 0.40696958949155611
Out[10]:
0.4069695894915561
In [11]:
draw_functions(-1,1,-1,1,100, [move_function(f,0,1,0,1)])
Out[11]:
In [12]:
test_function(f,0,1,0,1, exact_int, 30)
Out[12]:
In [13]:
test_exact_digits(f,0,1,0,1, exact_int, 30)
Out[13]:

Testy funkcji $\sin(x+y+1)$.

Możemy sprawdzić, że $\int_{-1}^1 \int_{-1}^1 \sin(x+y+1) dydx = 4sin^3(1)$.

In [14]:
f(x,y) = sin(x+y+1)
exact_int = 4*sin(1)^3
#draw_functions(-1,1,-1,1,100, [f])
Out[14]:
2.3832929463638224
In [15]:
#test_function(f,-1,1,-1,1,exact_int, 30)
In [16]:
#test_exact_digits(f,-1,1,-1,1,exact_int, 30)

Test funkcji $e^{x+y+1}$.

In [17]:
f(x,y) = exp(x+y+1)
exact_int = e*(e-1/e)^2
#draw_functions(-1,1,-1,1,100, [f])
Out[17]:
15.016852707441016
In [18]:
#test_function(f,-1,1,-1,1,exact_int, 30)
In [19]:
#test_exact_digits(f,-1,1,-1,1,exact_int, 30)

Testy funkcji $\sqrt{(x+1)(y+1)}$.

In [20]:
f(x,y) = sqrt((x+1)*(y+1))
exact_int = 32/9
draw_functions(-1,1,-1,1,100, [f])
Out[20]:
In [21]:
test_function(f,-1,1,-1,1,exact_int, 12)
Out[21]:
In [22]:
test_exact_digits(f,-1,1,-1,1,exact_int, 30)
Out[22]:

Test efektu Rungego.

In [23]:
f(x,y) = 1/(1 + 25x^2 + 25y^2)
exact_int = 0.4361934108122879
test_function(f, -1, 1, -1, 1, exact_int, 15)
Out[23]: